Theoretical Analysis of Subthreshold Oscillatory Behaviors in Nonlinear Autonomous 

Systems 
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We have developed a linearization method to investigate the subthreshold oscillatory behaviors in 
nonlinear autonomous systems. By considering firstly the neuronal system as an example, we show 
that this theoretical approach can predict quantitatively the subthreshold oscillatory activities, 
including the damping coefficients and the oscillatory frequencies which are in good agreement 
with those observed in experiments. Then we generalize the linearization method to an arbitrary 
autonomous nonlinear system. The detailed extension of this theoretical approach is also presented 
and further discussed. 
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Oscillatory behaviors are one of the most important 
features for the nonhnear systems and have attracted 
much attention in recent years [1, 2, 3|. Oscillation phe- 
nomena have been discovered in many areas of the chem- 
ical and biological sciences, such as chemical oscillation 
0i @] J rhythmic gene expression and metabolism Q , and 
particularly neuronal oscillations 0, H, . These oscilla- 
tory behaviors may have some observable influences on 
the responses of the systems, and play significant and 
often crucial roles, for example in neuronal system the 
oscillation activities are found to be important for in- 
formation processing [lol 11 | and cognitive perception 



[121 ll3L ll4L Il5l | . It is therefore of great interest to study 
the oscillatory behaviors in these nonlinear systems. Os- 
cillation theory has become an important part of con- 
temporary applied mathematics [3, [13] and has been 
well developed 18, [l^. However, the analysis of oscil- 
latory processes is still far from being complete. Up to 
date, little attention has been paid to the question how 
to systematically characterize the oscillatory behaviors, 
especially regarding the damping coefhcients and oscilla- 
tory frequencies, in these nonlinear systems. 

Generally, the dynamics of nonlinear systems are often 



described by some coupled differential equations [18|, 
While the stimuli to the system involve no variables that 
are differentiating (constant stimuli for example) we call 
the nonlinear system an autonomous one. In this Let- 
ter we attempt to approach the aforementioned issue by 
presenting a linearization method to theoretically inves- 
tigate the subthreshold oscillation behaviors in nonlinear 
autonomous systems. Firstly we consider the oscillatory 
activities in the neuronal system under subthreshold con- 
stant stimulus as an example. The theoretical approach 
enable us to obtain analytical solutions for the damp- 
ing coefficients and frequencies of oscillatory behaviors 
of membrane voltage. In addition, we also generalize 
the theoretical approach to an arbitrary nonlinear au- 
tonomous system. The detailed implementation of this 
theoretical framework is also presented and discussed. 

Neuronal systems are highly nonlinear excitable sys- 
tems. To elucidate the subthreshold oscillatory behav- 



iors in autonomous neuronal system, we consider the 
Hindmarsh-Rose (HR) neural model, as an example, un- 
der a subthreshold constant stimulus. The dynamics of 
HR neuron system could be described by the following 
equations (20j: 



dx/dt = y — z — ax^ + bx"^ + I(t), 
dy/dt = c — dx^ — y, 
dz/dt = r[s[x — xq) — z], 



(1) 



where x denotes the membrane potential, y the recov- 
ery variable and z a slow adaptation variable, and the 
parameters a, b, c, d, r, s and xq are the same as in 
Ref . [20| . If the membrane potential x{t) exceeds its 
threshold value Xth — 0, the neuron will result in a spike 
S{t) = 9{x{t) — Xth), where 9{x) is the Heaviside func- 
tion. I{t) represents the stimulus bias to the neuron, 
I{t) = Iq < Ic with Ic being the critical value for con- 
stant stimulus. 

It is well known that the HR neuron will be active with 
limit-circle firings [2l'| for a suprathreshold constant in- 
put, Iq > Ic in Eqs. (HI). If Iq < Ic, the system may un- 
dergo a temporal process of damping oscillation to reach 
a quiescent stable state, with the membrane potential 
x(t), the recovery variable y{t) and the slow adaptation 
variable z{t) being the stable value Xg, ys, Zs, respec- 
tively. These stable values can be obtained from Eqs. ([T]) 
by simply letting dx/dt = dy/dt = dz/dt — 0, yielding 



axl + (d — b)xl + sxs ~ {sxq + c) — Iq 



(2) 



and ys — c — dx'^, Zs — s{xs — xq). It is easy to prove 
that Eq. ([2]), with respect to Xs, has only one real root 
for any value of Iq, 

Xs - (-Q/2+VA)i/3 + (-Q/2-\^)i/^-(d-6)/3a (3) 



where A = {Q/2f 
2{d-bf/27a'^-s{d^ 
always write 



- [s/3a - 
6)/3a2 



{d - bf/Qa^f and Q = 
(sxo + c + /o)/a. One can 



x{t) = Xs + x{t), 

y{t) = ys + y{t), 
z{t)^Zs + S{t), 



(4) 
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Since /q < IcO, the variables x{t), y{t) and z{t) are small 
(approximately zero) for sufficiently large t; so one can 
linearize Eqs. H]), by utilizing Eqs. ^ into the following 
equation 



dx/dt 




'2bxs-3ax'^^ 1 -1" 




X 




dy/dt 




-2dxs -1 




y 


(5) 


dz/dt 




rs — r 




z 





The eigenvalue A of the matrix in Eqs. ^ yields 

{\+r)[{\+l){Wiaxl-2bxs) + 2dx,]+rs{\+l) =0 (6) 

When < /q < /c, Eq. ([6]) has one real root, which re- 
mains negative, and two conjugated complex roots, say 
\r ± iXi. These roots, as well as x^, depend on /q and 
can be obtained analytically through Eqs. ^ and ([6|): 

A,(/o) = -[(-e/2 + + (-e/2 - V5Y/^]/2 - a/3, 

A,(/o) = V3[(-e/2 - - (-e/2 + V5)i/3]/2, 

(7) 

where (5 = (e/2)2-K[(-aV3+/3)/3]3, e = 2aV27-a/3/3+ 
7, a = 3aa;^-2foa;s + r+l, /3 = 3a(l-hr)a;2-K2(cf-6)2;s-h 
r{s + l) and 7 = r[3aa;^ + 2((i— 6)a;s-|-s] with Xs obtained 
from Eq. 

The analytical results for the dependence of upon 
Jo are shown in figure 1. When /q increases from to 
the threshold value Jc, the negative A^ increases mono- 
tonically up to the critical value corresponding to the 
limit-circle firings behaviors observed in simulation 21 1. 
Be aware that there still exists damped behaviors even 
when the input to the neuron is a inhibitory one. The 
imaginary part A^ 7^ means that the membrane poten- 
tial x{t), as a function of the real time t, behaves a damp- 
ing oscillation with an intrinsic frequency fi = A^. The 
calculation results of fi versus /q are shown in figure 2. 
One can see that the intrinsic frequency fi of the damping 
oscillation is about 10 ~ 40 Hz, depending on Jq. It is no- 
table that the damping oscillation does exist even when 
Jo = 0, indicating the inherence of the oscillation. Obvi- 
ously, this frequency range is the same as the subthresh- 
old activity observed in experiments 0,122, H, 24, 25|. 

In fact, our theoretical analysis of subthreshold oscilla- 
tion activities in neuronal systems are naturally indepen- 
dent of the specific neuronal model used. The intrinsic 
oscillatory behaviors in other excitable neuronal mod- 
els, such as the Bonhoeffer van der Pol (BvP) model, 
FitzHugh-Nagumo (FHN) model and Hodgkin-Huxley - 
type (HH) one, can also be obtained analytically or nu- 
merically provided that the stimuli to these nonlinear 
systems are subthreshold constant ones. In the follow- 
ing we will describe a generic version of our theoreti- 
cal approach, for the case of an arbitrary nonlinear au- 
tonomous system under subthreshold stimuli, to investi- 
gate the intrinsic oscillatory behaviors with an emphasis 
on the damping coefficients and oscillatory frequencies. 

The dynamics of nonlinear autonomous systems are 
generally described by the differential equations of the 



first order or higher order for complex systems, and the 
latter can easily be reduced to several coupling differen- 
tial equations of first order. In this study, we suppose 
that an autonomous nonlinear system is governed by the 
following n coupling equations: 

dxi/dt = fi{xi,X2, ■■■Xn) + A, {i 1, 2, ...n) (8) 

where Ai are the constant stimuli to the system and make 
the nonlinear system an autonomous one. We simply let 
dxi /dt — and easily get the following n equations with 
real coefficients: 



fi{xi,X2, ■■■Xn) ^ 0, (i = l,2,...n) 



(9) 



In principle, there are up to n real roots for the above 
equations, indicating that the nonlinear autonomous sys- 
tem has n steady states at most. Generally, the Eqs. ([9|) 
may have both real roots and complex roots, the former 
of which implies the number of steady states while the 
latter should appear in pairs of conjugated complex due 
to the real coefficients of Eqs. (|9]). If all the roots of the 
Eqs. ^ are complex ones, we can conclude that the dy- 
namic of the nonlinear autonomous system will behavior 
divergently or chaotically and the nonlinear system is no 
longer convergent to any steady states, which is beyond 
our consideration in this study. 

For simplicity let us pay close attention to one of the 
steady state of system, Xs, where each variables of the 
system reaches their own steady value, x^s. During the 
convergent process of the system one can always write 



Xi — Xis 4- Xi , (2 — 1,2,. ■■Tl) 



(10) 



whereas small enough for sufficiently long time t. 

Inserting Eq. (fTO| into Eqs. ([8]) we can get 

dxjdt = /i(xis -I- Xi,X2s + X2, ■■■Xns + in) (« = 1, 2, ...n) 

(11) 

where the expression terms on the right side of Eqs. (fTTl 
can be further expanded at the steady state, Xs, up to 
the accuracy of first order Taylor expansion since Xi are 
sufficient small variables, into the following as: 

fi{xis + Xi,X2s + i2, ■■■Xns + in) 



Xi 



~ dfi 
^ ax2 



g^U^, (j = l,2,...n) 
(12) 

After inserting Eqs. ^ and Eqs. US]) into Eqs. HI]) we 
obtain n new differential equations expressed below in 
the form of the matrix. 



dxi/dt 




' 9/i 1 


dh\ 

dxi l^s 


dxi 1 dt 




dxi ^■^s 


dh\ 

dxi 


dxn/dt 






df^\ 

dxi l^s 



dfi I 
df^\ 

dXn 
dXn 



Xi 



(13) 



The eigenvalues of the nx n real matrix in Eqs. ([T3|) can 
thus be expressed as a function of steady values Xg . Since 
the coefficients of the matrix in Eqs. are real, the 
coefficients of the eigenfunction should also be real num- 
bers. In general the eigenvalues of the matrix with real 
coefficients should be either real numbers or in pairs of 
conjugated complex. For the case of real numbers, the 
eigenvalues would remain negative, corresponding to the 
convergence of the nonlinear system under subthreshold 
stimuli. Otherwise the system may undergo divergent 
or chaotic dynamics; For the case of pairs of conjugated 
complex, the real part of eigenvalue must be negative 
suggesting the damping factors of the system, while the 
imaginary part implies the frequencies of the oscillatory 
activities. As a rule, there are up to [^] components of 
oscillatory frequencies in this nonlinear autonomous sys- 
tem ([] means the biggest integer that no greater than 
what's in the square bracket). 

At the any other steady state beyond Xs , the same lin- 
earizaiton method could also be applied. The key of our 
theoretical approach is that the stimuli to the nonlinear 
autonomous system must be subthreshold ones, to keep 
the nonlinear system from divergence or chaos, so that 
the linearization method is always available. 

In summary, we have developed a linearization method 
to study the subthreshold oscillatory behaviors in non- 
linear autonomous systems. This theoretical approach 
enable us to predict the dynamics of subthreshold behav- 
iors and quantitatively describe the damping coefficients 
and frequencies of oscillatory behaviors, if exist. Given 
the case of nonlinear neuronal system under subthreshold 
constant stimuli, our theoretical analysis can give analyt- 
ical solutions about its oscillation behaviors that are con- 
sistent with experimental observations. Since oscillatory 
activities are surprisingly ubiquitous phenomena in a va- 
riety of nonlinear systems, the linearization method pro- 
vide us a possible theoretical approach to characterize the 
oscillatory behaviors in those autonomous systems. Re- 
garding the oscillatory behaviors in the non- autonomous 
nonlinear system, one have to seek other methods or so- 
lutions which are our interest of further study. 
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FIG. 1: Theoretical prediction of damping coefficients versus 
constant stimuli Jo. Ic is estimated rmmerically as the critical 
value of stimulus that make the membrane potential exceed 
its threshold value xth = after a sufficient long time, and 
chosen as 1.32. 
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FIG. 2: Theoretical prediction of oscillatory frequencies ver- 
sus constant stimuli Iq. Ic is chosen as the same as that in 
figure 1. 



